version 19
drop _all

* change the root folder according to your computer directory
global root = "D:\WorkingPaper-Series\Township_Light\Replication_Files"  

global dofiles = "$root\dofile"     
global working_data = "$root\workdata"
global tables = "$root\table"
global figures = "$root\figure"

cd "$working_data"


******************
*****Pre Set******
******************

global geo_ctrl = "topography_t* slope_t* roughness_t*"
global wth_ctrl = "gaez_t* precipitation_t* sunlight_t*"
global dist_ctrl = "county_dist_t* dist_hz_t* dist_coast_t*"

global full_ctrl = "$geo_ctrl $wth_ctrl $dist_ctrl"

global geo_ctrl_1 = "topography_t_1 slope_t_1 roughness_t_1"
global wth_ctrl_1 = "gaez_t_1 precipitation_t_1 sunlight_t_1"
global dist_ctrl_1 = "county_dist_t_1 dist_hz_t_1 dist_coast_t_1"

global full_ctrl_1 = "$geo_ctrl_1 $wth_ctrl_1 $dist_ctrl_1"

global geo_ctrl_2 = "topography_t_2 slope_t_2 roughness_t_2"
global wth_ctrl_2 = "gaez_t_2 precipitation_t_2 sunlight_t_2"
global dist_ctrl_2 = "county_dist_t_2 dist_hz_t_2 dist_coast_t_2"

global full_ctrl_2 = "$geo_ctrl_2 $wth_ctrl_2 $dist_ctrl_2"

global matchvar = "county_dist topography slope gaez sunlight dist_hz dist_coast"


global pre_ctrl = "light_t road_t pop_t"


********************
***** Figure 7 *****
********************

*****All Sample*****
use "township_panel_final.dta", clear

    forvalues k = 6(-1)1 {
       gen g_`k' = (period == -`k')
    }
    
    forvalues k = 0/5 {
         gen g`k' = (period == `k')
    }

replace g_1 = 0 

reghdfe road_ipol g_6-g4 treat_trend $full_ctrl, a(county_year idcode) cl(countycode)

*** Joint Significance of g_4-g_3: 0.117(0.056)
lincom g_4+g_3+g_2
*** Joint Significance of g0-g4: -0.094(0.057)
lincom g0+g1+g2+g3+g4

forvalues i = 1/11 {
    local m_`i' = e(b)[1,`i']
}

mat input mat1_fe = (`m_1',`m_2', `m_3',`m_4',`m_5',`m_6',`m_7', `m_8', `m_9', `m_10',`m_11')

mata st_matrix("V",sqrt(diagonal(st_matrix("e(V)"))))

mat V = V'

forvalues i = 1/11 {
    local m_`i' = V[1,`i']
}

mat input mat2_fe = (`m_1',`m_2', `m_3',`m_4',`m_5',`m_6',`m_7', `m_8', `m_9', `m_10',`m_11')


mat twfe = [mat1_fe\mat2_fe]
mat twfe = twfe'
svmat twfe

keep if twfe1 != .
keep twfe1 twfe2
rename (twfe1 twfe2) (coef std_err)

export delimited using "twfe_est_road_full.csv", replace

*****Sub-sample only include cohort >2016
use "township_panel_final.dta", clear

    forvalues k = 6(-1)1 {
       gen g_`k' = (period == -`k')
    }

    forvalues k = 0/5 {
         gen g`k' = (period == `k')
    }

replace g_1 = 0 

reghdfe road_ipol g_6-g4 treat_trend $full_ctrl if (first_year>=2016|first_year==0)&year>=2012, a(county_year idcode) cl(countycode)

*** Joint Significance of g_4-g_3: 0.092(0.047)
lincom g_4+g_3+g_2

forvalues i = 1/11 {
    local m_`i' = e(b)[1,`i']
}

mat input mat1_fe = (`m_1',`m_2', `m_3',`m_4',`m_5',`m_6',`m_7', `m_8', `m_9', `m_10',`m_11')

mata st_matrix("V",sqrt(diagonal(st_matrix("e(V)"))))

mat V = V'

forvalues i = 1/11 {
    local m_`i' = V[1,`i']
}

mat input mat2_fe = (`m_1',`m_2', `m_3',`m_4',`m_5',`m_6',`m_7', `m_8', `m_9', `m_10',`m_11')


mat twfe = [mat1_fe\mat2_fe]
mat twfe = twfe'
svmat twfe

keep if twfe1 != .
keep twfe1 twfe2
rename (twfe1 twfe2) (coef std_err)

export delimited using "twfe_est_road16.csv", replace
